Strong Coupling Solver for the Quantum Impurity Model 
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We propose a fast impurity solver for the general quantum impurity model based on the pertur- 
bation theory around the atomic limit, which can be used in combination with the local density 
approximation (LDA) and the dynamical mean field theory (DMFT). We benchmark the solver in the 
two band Hubbard model within DMFT against quantum Monte Carlo (QMC) and numerical renor- 
malization group (NRG) results. We find that the solver works very well in the paramagnetic Mott 
insulator phase. We also apply this impurity solver to the DMFT study of the anti-ferromagnetic 
phase transition in the unfrustrated Bethe lattice. The Neel temperature obtained by the fast im- 
purity solver agrees very well with the QMC results in the large Hubbard U limit. The method is a 
promising tool to be used in combination with the LDA+DMFT to study Mott insulators starting 
from first principles. 

PACS numbers; 71.10Fd,71.20.Be 



I. INTRODUCTION 



Recently, much effort has been devoted to develop 
methods for ab initio investigations of real materials with 
strongly correlated electrons. A most promising tool 
was build by combining the conventional first principle 
methods, such as the density functional theory in the 
local density approximation (LDA), with the newly de- 
veloped dynamical mean field theory (DMFT)[3|. Many 
numerical schemes such as LDA-hDMFT[|, LDA-h-hH, 
GW-|-DMFtQ has been proposed and applied to various 
systems 1^0. 

The application of DMFT to the real material re- 
quires a fast scheme to solve the generalized Anderson 
impurity model. Most of the impurity solvers, such as 
the iterative perturbation theory (IPT) 0, 01 1 the non- 
crossing approximation (NC A) j|| , the slave boson mean- 
field[lO| , the equation of motion method[T]|. exact diago- 
nalization_based methods |3 and quantum Monte Carlo 
(QMC) 13], have been developed for simplified multi- 
orbital Anderson models, usually assuming SU(N) sym- 
metry. While very few tools are available for the study 
of general Anderson impurity models generated by real- 
istic DMFT calculations. Therefore it is important to 
develop impurity solvers which can be used for a very 
general case. In the weak coupling limit when the sys- 
tem is in the metallic phase, the fluctuation exchange 
approximation(FLEX) jlJI has been proposed and imple- 
mented. On the other hand for the Mott insulator phase, 
when the local interaction term is very strong, we are still 
lack of such a general impurity solver which can treat 
models without SU(N) symmetry and containing general 
crystal fields and multiplet terms. Recent studies on the 



Mott insulators [1^, i.e. LaMnO^, V2O3, LaTiO^ and 
YTiOs 16], discovered a variety of phenomena includ- 
ing orbital order-disorder transition, charge density wave 
and anti-ferromagnetism. Therefore the first principle 
study on the Mott insulator material with or without long 
ranger order becomes a very important issue for both the 
theoretical understanding of these materials and material 
designing of these type of compounds. 



II. DERIVATION OF THE METHOD 

In this paper, we propose an impurity solver which 
is based on the perturbation theory around the atomic 
limit. We will consider the most general Anderson im- 
purity model generated by the LDA-I-DMFT calculation 
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where Hiocai is a very general atomic Hamiltonian and 
index a denotes the spin and orbital degree of freedom. 
Further, Hband stands for the conducting band which 
plays the role of a fermionic bath in DMFT calculations. 
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Our first step is to diagonalize the atomic Hamiltonian 
Hiocai by the exact diagonalization 



Hiocai = } E„,,\m){m\ 



(5) 



which can be done for any atom on modern computers. 
The hybridization term then takes the form 

Hv=Y^ VkapiF'^^Um' \m){m'\ckp + h.c. (6) 

kaf3 

where {F°')mm' — {'m-\fa\n^') , are the matrix elements 
of the operator in the local eigenbase. The atomic 
Green's function can be expressed by the eigenstates in 
the following way 
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where X,„ = e^^^^ jZ is the probability for the atomic 
state |m). The atomic self-energy can then be obtained 
by the inversion of a matrix 

= (iw + At - - G^''"'"'\iujf\ (8) 

This is the zeroth order self-energy in the expansion 
around the atomic limit. If we are able to compute the 
expansion of the Green's function in powers of the hy- 
bridization, i.e., G = G^"*""") + G(2) + 0{V^), we could 
also express the correction to the self-energy using the 
Dyson equation 

G = {iuj + fi-i-T.- Ay^. (9) 

To the lowest order in hybridization, we have 

^{atom) ^ ^(atom)^^ ^(2) ^{atom)^^ ^ -|- 0(X^^) 

(10) 



The expansion of the Green's function in the hybridiza- 
tion can be done by linked cluster expansion method fol- 
lowing Metzner et al. or using the auxiliary particle 
method T§\ or by straightforward expansion of the fol- 
lowing functional integral 



^ , , /^[/t/]4(0)/o(r)exp(-5) 



jD[fy]eM-S) 



where 



local 



a/3 



13 r0 



Jo 



dr / dr7l(T)A„^(r-r')/Mr')- 



(12) 

Expanding to the lowest order in A, one obtains two first 
order terms: a simple one body term from expanding the 
denominator and the complicated two-body term from 
expanding the nominator. The correction to the atomic 
Green's function takes the form 



G'fhiu) = G%'"^\iu)P Tr[G("*°™)A] + G\p{iLo) (13) 



where the two particle Green's function G^ is 



rP j-P pP 
Glyiu) = dT dTi drse'"" x (14) 
Jo Jo Jo 

Y.{TrUT)yyo)flin)fsiT,))o A^sin - r^). 



and the average (...)o is the atomic average, i.e., 
Tr{exp{— f3Hioc) ■ ■ ■ )/Z. Inserting the representation of 
the electron operator /„ = J2mm' ^mm'\'^)i'^'\ ^he 
above atomic average, one arrives at 
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The triple time integral can be done numerically or what lengthy algebra, the following expression for the 
analytically. In the latter case, one obtains after some- two-particle Green's function in the atomic limit 

J 
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where we used the notation Ei. 
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functions TZ-ys and Q-ys are 
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with the atomic probabilities given by = e ^^^/Z. 
On the real axis these two functions take the form 
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The constant in the first term of Eq. (|13|l can also be 
expressed by the above defined functions 

Tr[G("*°'»)A] = Y.iP'^oi{F-'^)w'Jl,s{Ei,Eo). (21) 

0,1,57 



n,s{Ei,E2) ^ Re{Xi A-,{Ei2)-X2 A+ (i?i2)} (18) C®' ED, 03, and O into Eq. m, 

i i 7'5^ ^ 7*^ ^^^J ^ ^ we finally obtain the self energy for the impurity mod- 

Qjs{ii^, El, E2) = Xi A^g{iuj + E12) + X2 A^g{iijj + E12) els requiring the impurity levels £, the interaction matrix 
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U and hybridization function A as an input. Thus we 
find a way to calculate the exact second order pertur- 
bation in the hybridization term. Implemented with the 
DMFT self consistent condition described in reference 
0, this method can be used as a very efhcient impurity 
solver in the DMFT study of the multi-orbital systems 
with very complicated local interactions. 



III. BENCHMARK 

To test this impurity solver, we calculated the Green's 
function for the two band Hubbard model at half fill- 
ing and compared with the results obtained by the QMC 
solver. We choose the temperature to be 0.125, where the 
QMC result is quite reliable. We use the semi-circle den- 
sity of states (DOS) and set the half bandwidth D = 1 
as the unit of energy. The metal-insulator transition has 
been determined by QMC at Uc = 3.5. First let's com- 
pare the results for U — 6, where the system is on the 
insulator side. We found that for frequency higher than 
u! = which is the highest energy scale in the prob- 
lem, all the results obtained from three different scheme 
(QMC, 0-th order atomic expansion, in which we sim- 
ply use the atomic self energy defined by Eq. ((HJ and 
the present solver) fall onto a single curve, which indi- 
cates that both the 0-th order atomic expansion and the 
present solver can capture the correct high energy limit. 
While for the frequency lower than U, the result obtained 
by 0-th order atomic scheme show clear deviations from 
the QMC data, including the the second order strong cou- 
pling perturbation correction implemented by the present 
solver gives excellent results as shown in Fig.^ The situ- 
ation is similar for U = A which is close to the Mott tran- 
sition point but still on the insulator side. In this case, 
the deviation between the 0-th order result and the QMC 
result becomes quite large at low frequency, as shown in 
Fig. 121 For example, the relative error reaches 53% at the 
first Matsubara frequency. Again the error is corrected 
by turning on the second order perturbation around the 
atomic limit in the present solver. Based on the above 
comparison, we can draw the conclusion that the second 
order perturbation in the hybridization term works very 
well in the Mott insulator phase. 

Another good test for this impurity solver is the anti- 
fcrromagnetic(AF) order in the half filled single band 
Hubbard model on a bipartite lattice. In the large U 
limit, it is well known that in such a case up to the sec- 
ond order perturbation in t/U, the Hubbard model can 
be mapped to an anti-ferromagnetic Heisenberg model 
with exchange energy J cx jU . Since the exchange en- 
ergy is the only energy scale in the problem, the Neel 
temperature must also be proportional to 1/U in large 
U limit. Although the physics behind is quite straight- 
forward, it is not reproduced by EOM or IPT methods 
within the framework of DMFT. Although the Hartree- 
Fock approximation can also obtain the correct AF order 
in the ground state, it predicts the wrong Neel temper- 
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FIG. 1: The comparison of QMC and DMFT with atomic 
solver for two band Hubbard model with U = 6, n — ,T = 
0.125 and half band width D = 1. 



-10 
-12 



r 






— present solver 






- - ■ - - 0th order atomic 






-•-DMFT+QMC 











10 12 



FIG. 2: The comparison of QMC and DMFT with atomic 
solver for two band Hubbard model with U = 4, jj. — ,T = 
0.125 and half band width D = 1. 



ature to be of order U instead of J, as shown in Fig. O 
Since the impurity solver we proposed here can include 
exactly the second order correction to the atomic limit, 
we expect that it can reproduce the correct Neel tem- 
perature in the large U limit. For this purpose, we 
solved the single band Hubbard model at half filling on 
the unfrustrated Bethe lattice and compared the results 
with the Hartree-Fock approximation, DMFT-I-IPT and 
DMFT+QMC in FigEl Since the self energy obtained 
by the present impurity solver vanishes when U — 0, 
we can also get the correct result for the non-interacting 
case. That is why the Neel temperature obtained by 
the present solver goes down in the small U limit. A 
very good agreement between our results and the QMC 
results is found for U/D > Uc, where Uc = 3.5 is ap- 
proximately the Uc2 in the Mott transition in the param- 
agnetic phase We also show the results obtained by 
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FIG. 3: The comparison of the Neel Temperature obtained 
by QMC,IPT, Hatree-Fock and atomic solver for half filled 
single band Hubbard model with half band width D = 1. 

iterative perturbation theory and Hartree-Fock approx- 
imation in figure O which are far away from the QMC 
results for almost the whole range of U. The remark- 
ably good agreement between our results and the QMC 
data in the large U limit indicates that combined with 
LDA, this simple impurity solver can be used to carry 
out the first principle calculation of the ordering tem- 
perature of the materials which have spin/orbital long 
range order in the ground states in the framework of 
LDA+DMFT. Compared with the model Hamiltonian 
studies, in which the super exchange processes are con- 
sidered by the Heisenberg model, the LDA-I-DMFT ap- 
proach has two advantages. The first one is that unlike 
the Heisenberg model, which can only capture the low 
energy physics, the LDA-fDMFT can capture not only 
the low energy physics like the long range spin/orbital 
order but also the high energy physics like the Hubbard 
bands. Besides that, since the effective bath in DMFT 
does not come only from the nearest-neighbor sites, the 
LDA+DMFT calculation can include the long range cou- 
pling between the local spins in a natural way. 

For the paramagnetic phase of an SU(N) Anderson im- 
purity model, the self-energy of the present solver is sim- 
ple enough to be written as a closet expression. It be- 
comes particularly simple in the case of half-filing where 
it takes the form 

To get the DMFT solution, we also need the DMFT 
self-consistency condition. For the Bethe lattice, it is sim- 
ply given by A = t'^G therefore the DMFT local Green's 
function takes the form 



Gdmft{z) = - s^x^-AT{zy) (23) 



FIG. 4: The spectral function of the one band Hubbard model 
on the Bethe lattice for various values of U compared with the 
NRG results (taken from Ref. 

where t(z) = + x - ^^^^ ^ ^ 

sign(Im(a;^ — 4r(2:)^)). 

The spectral function that corresponds to Eq. H23() is 
plotted in Fig. 01 for various values of J7, ranging from 
U = to J7 = 6. For comparison, we also displayed the 
NRG results at U = A. It is clear, that the present solver 
misses the Kondo effect and therefore should be used only 
in the insulating state, i.e., when the spectral function 
has a finite gap. At U = 4 we can see that the width 
of the gap as well as the position of the Hubbard bands 
is very close to the NRG results. Note, however, that 
the NRG has a finite resolution at high frequencies and 
therefore does not provide very precise Hubbard bands 
either. They usually tend to be slightly rounded. 

It is well known that deep in the insulating state, the 
width of the Hubbard bands as well as the shape of the 
Hubbard bands has to be the same as the non-interacting 
density of states. As one can see in Fig. 0] the width is 
indeed 2D and they become more and more semicircular 
as U increases. Note that the zeroth order approxima- 
tion (expressed by Eq. (jHJ) gives for factor of 2 narrower 
Hubbard bands. Finally, the critical U, at which the gap 
closes, is 2^/3 and is reasonable close to the exact upper 
critical U being Uc2 2.97. 

It should be noted that the perturbation theory in the 
hybridization is a singular perturbation for the metallic 
system, therefore any finite order perturbation can not 
offer a qualitatively proper description of the system. If 
one anyway applies the present solver to the metallic sys- 
tem, the local Green's function develops a pole in the 
complex plane which does not lie on the real axis. This 
pole indicates a tendency toward the formation of a sin- 
gularity at the Fermi level. Namely, some weight is miss- 
ing under the Hubbard bands and the spectral function 
develops a V shaped cusp at the Fermi level. 

To avoid the causality problem in the metallic state, 
one might rewrite the self-energy in the continued frac- 
tion representation which has the same lowest order term 
in expansion in A. For the half- filled one band model. 
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FIG. 5: The DMFT spectral function obtained by using self- 
energy from Eq. 1241 compared with the NRG results. 

the following self-energy can be constructed in this way 
E(z) ^ (-\ 1-—. (24) 

The DMFT spectral function that corresponds to this 
self-energy is plotted in Fig. [S] In this approximation, 
the system is metallic for U < V3, however, the metal- 
lic state is not Fermi-liquid and therefore the spectral 
function does not reach the unitary limit at = 0. How- 
ever, the causality problem is avoided and the impurity 
solver does not break down in the metallic state. The 
agreement in position and shape of the Hubbard bands 
is also improved in the insulating state of the system. 
Note almost perfect agreement between NRG and the 
present solver spectral function at [/ = 4. As an effi- 
cient impurity solver in the strong coupling limit, the 
present solver should be compared with two other im- 
purity solvers, which are commonly used in this limit, 
namely the equation of motion methodfEOMl [lll| and 
NCA 9]. Unlike the method proposed in this paper, the 
EOM method requires a self consistent procedure to ob- 
tain the Green's function on the impurity site. A gen- 
eral impurity model generated by LDA-fDMFT usually 
contains a large numbers of orbits and very complicated 



non-SU(N) like local interaction. Therefore in the EOM 
method one has to solve the self consistent equations with 
very large number of parameters, which is almost impos- 
sible numerically when the orbital number reaches 14 like 
in systems with one open / shell per unit-cell. The NCA 
method suffered from the same problem when the orbital 
number becomes large and the system is away from the 
SU(N) symmetry since the number of atomic states, and 
therefore pseudo-particles, grows exponentially. Com- 
pared with other impurity solvers, the main advantage 
of the present solver is that no self consistent loop is 
required to solve the impurity problem and the local in- 
teraction can be very general including local Coulomb re- 
pulsion, Hunds coupling, spin-orbital coupling and pair 
hopping term. The computational time of the present 
impurity solver still grows as where N is the num- 
ber of atomic states, since one needs to carry out four 
sums in Eq. H16|l . Note, however, that the matrix 
has many zero elements and if one takes into account the 
conservation of spin and particle number, the computa- 
tional time can be considerably reduced. At present, for 
a general / system, it takes only few seconds on modern 
computers. 



IV. CONCLUSIONS 

To summarize, this paper presents a new impurity 
solver based on second order perturbation theory in the 
hybridization around the atomic limit. The strength of 
the approach lies in its generality and its speed. It can be 
applied to systems with very complicated atomic config- 
urations, general Coulomb repulsion and spin-orbit cou- 
pling and does not require self-consistency therefore it 
can be efficiently applied to a systems with open d or f 
shells. It is however limited to integer filling and large 
U, which is the regime particular important to study 
transition-metal compounds in the Mott insulating state 
with or without orbital or magnetic long-range order. 

ACKNOWLEDGEMENT: This work was supported 
by the NSF under grant DMR-0096462. 



[1] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, 
Rev. Mod. Phys. 68, 13 (1996). 

[2] V. Anisimov, A. Poteryaev, M. Korotin, A. Anokhin, 
and G. Kothar, J. Phys.: Cond. Matt. 9, 7359 (1997); 
K. Held, I.A. Nekrasov, G. Keller, V. Eyert, N. Bliimer, 
A.K. McMahan, R.T. Scalettar, Th. Pruschke, V.I. Anisi- 
mov, and D. VoUhardt, Psi-k Newsletter #56 (April 
2003), p. 65; A. I. Lichtenstein, M. I. Katsnelson, and 
G. Kotliar, in Electron Correlations and Materials Prop- 
erties 2, ed. A. Gonis (Kluwer,NY) cond-mat/0211076 . 

[3] A.I. Lichtenstein and M. I. Katsnelson, Phys. Rev. B 57, 
6884 (1998). 

[4] S. Biermann, F. Aryasetiawan, and A. Georges, Phys. 



Rev. lett. 90, 086402 (2003). 
[5] X. Dai, S.Y. Savrasov, G. Kotliar, A. Mighori, H. Led- 
better, and E. Abrahams, Science 300, 953 (2003); S.Y. 
Savrasov, G. Kothar, and E. Abrahams, Nature 410, 793 
(2001). 

[6] K. Held, G. Keller, V. Eyert, D. VoUhardt, and V.I. 

Anisimov, Phys. Rev. Lett. 86, 5345 (2001). 
[7] M.J. Rozenberg, G. Kothar, and X.Y. Zhang, Phys. Rev. 

B 49, 10181 (1994). 
[8] S.Y. Savrasov, V. Oudovenko, K. Haule, D. Villani, and 

G. Kotliar Phys. Rev. B 71, 115117 (2005). 
[9] M. Jarrell, and Th. Pruschke, Phys. Rev. B 49, 1458 

(1994); K. Haule, V. Oudovenko, S.Y. Savrasov, and G. 



7 



Kotliar, Phys. Rev. Lett. 94, 036401 (2005). 
[10] G. Kotliar and A.E. Ruckenstein, Phys. Rev. Lett. 57, 
1362 (1986). 

[11] H. O. Jeschke and G. Kotliar Phys. Rev. B 71, 085103 

(2005); Jian-Xin Zhu, R. C. Albers, J. M. Wills, 

cond-mat/0409215 
[12] M. Capone, M. Fabrizio, C. Castellani, and E. Tosatti, 

Phys. Rev. Lett. 93, 047001 (2004); A. Koga, N. 

Kawakami, T. M. Rice, and M. Sigrist, Phys. Rev. Lett. 

92, 216402 (2004). 
[13] L. Laloux, A. Georges, and W. Krauth, Phys. Rev. B 50, 

3092 (1994). 

[14] L. Chioncel, L. Vitos, LA. Abrikosov, J. KoUr, M. L 
Katsnelson, and A.I. Lichtenstein, Phys. Rev. B 67, 
235106 (2003); V. Drchal, V. Janis, J. Kudrnovsk, V.S. 
Oudovenko, X. Dai, K. Haule, and G. Kotliar, J. Phys.: 



Cond. Matt. 17, 61 (2005). 

[15] Y. Tokura, and N. Nagaosa, Science 288, 462 (2000). 

[16] E. Pavarini, S. Biermann, A. Poteryaev, A.I. Lichten- 
stein, A. Georges, and O.K. Andersen, Phys. Rev. Lett. 
92, 176403 (2004). 

[17] W. Metzner, Phys. Rev. B 43, 8549 (1991); T.D. 
Stanescu, and G. Kotliar, Phys. Rev. B 70, 205112 
(2004). 

[18] K. Haule, S. Kirchner, J. Kroha, and P. Wolfle, Phys. 

Rev. B 64, 155111 (2001). 
[19] A. Georges and W. Krauth, Phys. Rev. B 48, 7167 

(1993); M.J. Rozenberg, G. Kotliar, and X.Y. Zhang, 

Phys. Rev. B 49, 10181 (1994). 
[20] R. Bulla, .cond-mat/0412314^ 



